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ABSTRACT 

A new Navier-Stokes algorithm for use on unstructured triangular meshes is presented. Spatial 
discretization of the governing equations is achieved using a finite-element Galerkin approxi- 
mation, which can be shown to be equivalent to a finite-volume approximation for regular 
equilateral triangular meshes. Integration to steady-state is performed using a multi-stage time- 
stepping scheme, and convergence is accelerated by means of implicit residual smoothing and 
an unstructured multigrid algorithm. Directional scaling of the artificial dissipation and the 
implicit residual smoothing operator is achieved for unstructured meshes by considering local 
mesh stretching vectors at each point. The accuracy of the scheme for highly stretched triangu- 
lar meshes is validated by comparing computed flat-plate laminar boundary-layer results with 
the well known similarity solution, and by comparing laminar airfoil results with those 
obtained from various well-established structured quadrilateral-mesh codes. The convergence 
efficiency of the present method is also shown to be competitive with those demonstrated by 
structured quadrilateral-mesh algorithms. 
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Applications in Science and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23665. 
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1. INTRODUCTION 

The use of unstructured triangular meshes in two dimensions and tetrahedral meshes in 
three dimensions has proven valuable for computing inviscid compressible flow about complex 
geometries [1,2,3]. Unstructured meshes also provide a natural setting for the use of adaptive 
meshing, which has been shown to provide large increases in efficiency and accuracy [3,4], 
However, triangular and tetrahedral meshes have seldom been employed for computing viscous 
flows. Solutions of the full Navier-Stokes equations on triangular meshes can be found in the 
literature [3, 5, 6, 7]. However, these are often limited to low-Reynolds-number flows, and/or the 
accuracy and efficiency of these methods is inferior to that of existing quadrilateral mesh 
solvers. Consequently, numerous attempts at solving viscous flows for non-simple 
configurations have resorted to hybrid structured-unstructured meshing strategies, where struc- 
tured quadrilateral meshes are employed in the viscous regions, and unstructured meshes are 
employed in the inviscid regions. While such strategies have proven valuable for computing 
flows over various types of configurations [8], they lack the generality required for arbitrarily 
complex geometries. The use of completely unstructured meshes in both viscous and inviscid 
flow regions, as proposed in [9], holds the promise of producing a more general and flexible 
method for computing viscous flows over truly arbitrary configurations, while providing an 
ideal setting for the use of adaptive meshing techniques through the viscous layers, as well as 
in regions of strong viscous-inviscid interactions. 

For high-Reynolds-number flows over streamlined bodies, viscous effects are confined to 
thin boundary-layer and wake regions. As the Reynolds number increases, the viscous regions 
generally become thinner, and the gradients in the normal direction within these regions 
increase. To accurately resolve such flows, a small mesh spacing is required within the viscous 
regions. Since the flow gradients are predominantly in the normal direction, it proves economi- 
cal to refine the mesh only in this direction, retaining a large spacing in the tangential direc- 
tion. This approach is often employed for quadrilateral meshes, and may result in rectangular 
cells in the viscous regions with aspect ratios up to 10,000:1 for Reynolds numbers of 10 mil- 
lion. Clearly, for such cases, refinement in both normal and tangential directions would be 
prohibitively expensive. Thus, a directional refinement or stretching of the mesh must be 
employed for triangular meshes as well. This results in highly skewed triangles in the viscous 
regions, which may potentially degrade the accuracy and efficiency of the scheme. 

In this work, a Navier-Stokes solver for unstructured triangular meshes is described. A 
previously developed unstructured multigrid algorithm [10] is employed to accelerate the con- 
vergence of the solution to steady-state. Our objective is to demonstrate that by carefully tailor- 
ing the scheme for directionally stretched meshes, accurate and efficient solutions can be 
obtained which are competitive with those produced by current state-of-the-art structured 
quadrilateral-mesh Navier-Stokes solvers. While the solutions presented in this paper consist of 
laminar flow cases computed on regular stretched triangulations, this represents the necessary 
first step for validating the proposed algorithm, and establishes the feasibility of computing 
viscous flows on highly stretched unstructured meshes which contain a smooth variation of ele- 
ments through the viscous and inviscid regions. 

2. DISCRETIZATION OF THE GOVERNING EQUATIONS 

In conservative form, the full Navier-Stokes equations read 
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where w is the solution vector and f c and g c are the cartesian components of the convective 
fluxes 
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In the above equations, p represents the fluid density, u and v the x and y components of fluid 
velocity, E the total energy, and p is the pressure which can be calculated from the equation of 
state of a perfect gas 
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where cr represents the stress tensor, and q the heat flux vector, which are given by the consti- 
tutive equations for a Newtonian fluid 
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Y is the ratio of specific heats of the fluid, A/„ the freestream Mach number. Re. the Reynolds 
number based on the airfoil chord, and Pr the Prandti number. The coefficient of viscosity p 
varies with the temperature of the fluid, and is calculated as 


p = K T on (6) 

where K is a constant Equation (1) represents a set of partial differential equations which 
must be discretized in space in order to obtain a set of coupled ordinary differential equations, 
which can then be integrated in time to obtain the steady-state solution. 

The spatial discretization procedure begins by storing flow variables at the vertices of the 
triangles. The stress tensor o and the heat flux vector q must be calculated at the centers of the 
triangles. This is achieved by computing the required first differences in the flow variables 
(from equations (5)) at the triangle centers. For a piecewise linear approximation of the flow 
variables in space, the first differences are constant over each triangle, and may be computed 
as 
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where the summation over k refers to the three vertices of the triangle. The flux balance equa- 
tions are obtained by a Galerkin finite-element type formulation. The Navier-Stokes equations 
are first rewritten in vector notation 


V.F„ 


di Re. 

where the bold typeset denotes vector quantities, 
integrating over physical space yields 

dxdy + Jj\j>V.F c dxdy = ^ - fj^V.F, dxdy 

Integrating the flux integrals by parts, and neglecting boundary terms gives 

dxdy = J|F C .V<|> dxdy - jjF v .V<|> dxdy 


(9) 


Multiplying by a test function (|>, and 
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In order to evaluate the flux balance equations at a vertex P, <t> is taken as a piecewise linear 
function which has the value 1 at node P, and vanishes at all other vertices. Therefore, the 
integrals in the above equation are non-zero only over triangles which contain the vertex P, 
thus defining the domain of influence of node P, as shown in Figure 1. To evaluate the above 
integrals, we make use of the fact that and <J> y are constant over a triangle, and may be 
evaluated as per equations (7) and (8). The convective fluxes F e are taken as piecewise linear 
functions in space, and the viscous fluxes F v are piecewise constant over each triangle, since 
they are formed from first derivatives in the flow variables. Evaluating the flux integrals with 
these assumptions, one obtains 
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where the summations are over all triangles in the domain of influence, as shown in Figure 1. 
A M represents the directed (normal) edge length of the face of each triangle on the outer boun- 
dary of the domain, F* Ff are the convective fluxes at the two vertices at either end of this 
edge, and Ft is the viscous flux in triangle e. If the integral on the left hand side of equation 
(12) is evaluated in the same manner, the time derivatives become coupled in space. Since we 
are not interested in the time-accuracy of the scheme, but only in the final steady-state solution, 
we employ the concept of a lumped mass matrix. This is equivalent to assuming w to be con- 
stant over the domain of influence while integrating the left hand side. Hence, we obtain 
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where the factor of 1/3 is introduced by the integration of <t> over the domain, and Sl p represents 
the surface area of the domain of influence of P. For the convective fluxes, this procedure is 
equivalent to the vertex finite- volume formulation described in [1,10]. For the viscous fluxes, 
in the case of equilateral triangles, this formulation can be shown to be equivalent to a finite- 
volume formulation where the control volume is taken as the hexagonal cell formed by joining 
the centroids of all triangles with a vertex at P, as shown in Figure 2. For a smoothly varying 
regular triangulation, the above formulation is second-order accurate. 
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3. ARTIFICIAL DISSIPATION 

In principle, the physical viscous terms of the Navier-Stokes equations are capable of 
providing the numerical scheme with the dissipative property necessary for stability and captur- 
ing discontinuities. However, for high-Reynolds-number flows, this can only be achieved by 
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resorting to extremely small mesh spacings throughout the domain. Thus, in practice, it is 
necessary to introduce artificial dissipative terms to maintain stability in the essentially inviscid 
portions of the flow field, and to efficiently capture discontinuities. These additional dissipative 
terms must be carefully constructed to ensure that the accuracy of the scheme is preserved both 
in the inviscid region of the flow field where the convective terms dominate, as well as in the 
boundary-layer and wake regions where the artificial dissipation terms must be much smaller 
than the physical viscous terms. Previous Navier-Stokes solutions on highly stretched struc- 
tured meshes [11,12,13] have demonstrated the need for different scalings of the artificial dissi- 
pation terms in the streamwise and normal directions within the regions of viscous flow. How- 
ever, for unstructured meshes, directional scaling is significantly more difficult to achieve since 
no mesh coordinate lines exist. In fact, unstructured meshes have traditionally been considered 
to be truly multi-dimensional isotropic constructions with no preferred directions. However, as 
stated previously, the efficient solution of high-Reynolds-number viscous flows requires the use 
of meshes with highly stretched elements in the boundary-layer and wake regions, since these 
physical phenomena are highly directional in nature. For such meshes, even in the unstructured 
case, a direction and magnitude of stretching can be defined for each mesh point, as shown in 
Figure 3. This stretching vector, denoted as s, need not necessarily line up with any of the 
mesh edges. For the meshes employed in this paper, which are directly derived from structured 
quadrilateral meshes by splitting each quadrilateral into two triangles, the stretching magnitude 
and direction may be taken as the aspect-ratio and the major axis of the generating quadrila- 
teral element for each triangular element respectively. In the more general case, the generation 
of directionally stretched unstructured meshes [9,14] requires the definition of local stretching 
vectors throughout the flow field. These can in turn be used to scale the dissipation terms. It is 
important to note that these stretching vectors represent grid metrics which do not depend on 
the flow solution. 

The artificial dissipation operator on unstretched unstructured meshes has previously been 
constructed as a blend of an undivided Laplacian and biharmonic operator in the flow field. 
Since the biharmonic operator may be viewed as a Laplacian of a Laplacian, the dissipation 
operator may be reformulated as a global undivided Laplacian operating on a blend of the flow 
variables and their second differences : 

D(W) = flattta+Kjy] (14) 

where 

u = - iqV 2 * (15) 

In the above equations, fl represents the area of the control volume, which is of order Ax 2 , and 
V 2 w denotes the undivided Laplacian of w. The first term in the above equation constitutes a 

relatively strong first-order dissipation term which is necessary to prevent unphysical oscilla- 

tions in the vicinity of a shock. To preserve the second-order accuracy of the scheme, this term 
must be turned off in regions of smooth flow. This is accomplished by evaluating tc^ at mesh 
point i as 

El Pk ~ Pi 

(Kali = k 2 (16) 

2>* + Pi 
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Hence k'z is proportional to an undivided Laplacian of the pressure, which is constructed as a 
summation of the pressure differences along all edges meeting at node i, as depicted in Figure 
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3. This construction has the required property of being of order unity near a shock and small 
elsewhere. k 2 is an empirically determined coefficient which is taken as 0 for subcritical flows, 
and as 1/2 for transonic and supersonic flows. In equation (14), the overall scaling of the dissi- 
pation is performed by the factor a, which has previously been taken as proportional to the 
maximum eigenvalue of the Euler equations for inviscid flow calculations [4]. Directional 
scaling of the dissipation may thus be achieved by replacing equation (14) by 


D(w) = n [ctju^ + (17) 

where <Xi and ctj represent the different scalings in the $ and q directions respectively. Here, £ 
denotes the direction of the mesh stretching, and q the direction normal to £. Appropriate 
expressions for a : and ct 2 remain to be determined, as well as the discretization procedure for 
the above operator on unstructured meshes. 

On structured meshes, the dissipation is often scaled by the maximum eigenvalue of the 
Euler equations in each mesh coordinate direction, which is given by 

A*=[|«| + c]Aq X^Hvl + cJAS (18) 

where u, v, and A!;, Aq represent the local fluid velocity components and the mesh spacing 
respectively, in computational space, and c denotes the local speed of sound. However, for 
efficient multigrid convergence, a more even distribution of the dissipation is required in the 
two mesh coordinate directions, and the above scaling is replaced by [1 1,13]: 

h, = m h X n = <Kr- 1 )X n (19) 

where 


<Kr) = 1 + r 2 * 
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and 
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On unstructured meshes, we begin by constructing an isotropic value of the maximum eigen- 
value at each mesh point as 


X = ^|u x dl| + c|dl| (22) 

where the integration is performed around the boundary of the control volume for the particular 
mesh point being considered, and the bold typeset denotes vector quantities. The discrete 
approximation to the above integral yields the final form for X 


X = X I u AB&yAB ~ v ab^ x ab I + CxsVAx * \b + Ay 2 ab (23) 

where Ax^ and A represent the x and y increments along the outer edge AB of element e, as 
shown in Figure 1, and u^, v^, and c w represent averaged values along the edge AB. By 
considering the equivalent integration around the control volume on a structured quadrilateral 
mesh, it can be seen that X approximates the sum of the eigenvalues in the two space dimen- 
sions, i.e. 


X = V Xt, (24) 

Furthermore, the magnitude of the stretching vector s on the unstructured mesh can be con- 
sidered to be closely related to the cell aspect-ratio. Thus, by analogy with the structured mesh 
case 
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where s represents the magnitude of s , and the second approximation assumes that the magni- 
tude of the speed of sound c is much greater than the streamwise and normal velocities u and v 
in the viscous flow regions. Thus, equations (24) and (25) permit an estimate of the values of 
the maximum eigenvalues in the directions parallel and normal to the local mesh stretching 
vector, given the values of X and s. From equations (19) through (21), the cti and a 2 
coefficients of equation (17) are constructed as 

a, = m -^X 02 = 4 **- 1 ) ~x ( 26 ) 

Next, the discretization of the scaled Laplacian of equation (17) on unstructured meshes must 
be considered. Previously, for inviscid flows [1,4], the unsealed Laplacian of equation (14) was 
approximated as an accumulated edge difference in computational space, ie. 

“x* + Uyy = 77 £ [m* - «,] (27) 

“ *-i 

where k=l,...,n represents the n neighbors of node i, and the difference is taken along all edges 
meeting at node i. For a cartesian grid, this reduces to the familiar five point Laplacian finite- 
difference formula. Equation (17) can easily be approximated on a cartesian mesh aligned with 
the £, and q coordinate directions, simply by multiplying the constructed second differences in 
the \ and q directions by aj and a 2 respectively. For more general unstructured mesh topolo- 
gies, a finite- volume approximation to the scaled Laplacian of equation (17) yields the follow- 
ing discretization formula for the dissipation operator (see Appendix A) 


n 

D(Wi) = £2 [ + a 2 «^, ] = £[«*- UiHcqcosfy + a 2 sin 2 0J (28) 

*= i 

where 0* represents the angle between the kth mesh edge at node i, and the principal stretching 
direction as shown in Figure 3. From the above equation, it can be seen that if the kth mesh 
edge coincides with the £ or q directions, then the difference along that edge is multiplied by 
cti or a 2 respectively, and if cti = o t 2 , then the above discretization reduces to the isotropic 
accumulated edge difference previously employed. Since in practice ai and vary throughout 
the mesh, equation (28) is replaced by 
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Ak = a lt cos 2 0* + Oik sin 2 0 t (30) 

and the i and k subscripts refer to variables evaluated at nodes i and k, thus ensuring a conser- 
vative formulation of the dissipation operator. 


4. INTEGRATION TO STEADY-STATE 

The discretization of the spatial derivatives transforms equation (1) into the set of coupled 
ordinary differential equations 

dwi 

n ‘~dT + lQ(Wi) " D(m,, ' )] = °’ i=h2,3 n (31) 

where n is the number of mesh nodes. The residual Q(w) represents the discrete approximation 
to the convective fluxes. D(w) now represents the dissipative terms, i. e. the discrete approxi- 
mation to the viscous fluxes, as well as the artificial dissipation terms. These equations are 
integrated in time using a five-stage hybrid time-stepping scheme given by 

w (0) = w" 



( 32 ) 
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w" represents the value of the solution vector at the nth time step, and wW represents the value 
at the qth stage within a time step. The dissipative operator D(w) is evaluated only at the first, 
third, and fifth stages of the scheme, and is employed to construct the subscripted D q operator 
which represents a linear combination of present and previous evaluations of D(w). This 
scheme represents a particular case of a large class of multi-stage time-stepping schemes where 
the coefficients are chosen in order to maintain good stability properties when the viscous 
terms are dominant, and to ensure large damping of high-frequency errors, which is crucial for 
a rapidly convergent multigrid method [11]. The values of these coefficients are taken as 


and 


P = 0.56 y= 0.44 


cti = 1/4 (*2 = 1/6 a 3 = 3/8 04 = 1/2 05 = 1 


4.1. Local Time-Stepping 

Convergence to the steady-state solution may be accelerated by sacrificing the time accu- 
racy of the scheme, and advancing the equations at each mesh point in time by the maximum 
permissible time step in that region, as determined by local stability analysis. Stability limita- 
tions due to both the convective and diffusive characters of the Navier-Stokes equations must 
be considered. The local time step is thus taken as 


A/ = CFL ( 


to e + A/, 


■) 


(33) 


where CFL is the Courant number for the particular time-stepping scheme, and A t e and Af v 
represent the individual convective and viscous time-step limits respectively. The convective 
time-step limit has previously been derived for Euler solutions on unstructured meshes [4], and 
is given by 


At ' = X < 34 > 

where X e , previously denoted simply as X, represents the maximum eigenvalue of the inviscid 
equations averaged around the boundary of the control volume, as given in equation (23), and 
fl denotes the area of the control volume. The viscous time-step limit is taken as 


Ar y 



(35) 


- 8 - 


where K v is an empirically determined coefficient which determines the relative importance of 
the viscous and inviscid time-step limits in the final expression, and has been taken as 0.25 in 
this work. K represents the maximum eigenvalue of the diffusive operator of the Navier-Stokes 
equations, averaged about the boundary of the control volume. For the structured mesh case, 
A^ and X*, in the two mesh coordinate directions have been derived in [11]. For example, 




Re Q 


with a similar expression for 
may thus be approximated as 


Pr P 
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If the cross terms are neglected, K for unstructured meshes 


X. — 


Re Q 
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where the integration is performed along the boundary of the control volume. In discrete form, 
the expression for A, becomes 


where and p^ represent averaged values of viscosity and density along the outer edge AB 
of each element e. (c.f. Figure 1) 


4.2. Implicit Residual Smoothing 

The stability range of the basic time-stepping scheme can be increased by implicitly 
smoothing the residuals. Thus the original residuals R may be replaced by the smoothed residu- 
als R^by solving the implicit equations 

Ri=Ri + zV% (39) 

at each mesh point i, where V 2 /?, represents the undivided Laplacian of the residuals, and e is 
the smoothing coefficient For highly stretched structured meshes, the use of individual 
smoothing coefficients in the £ and q mesh coordinate directions which vary locally throughout 
the mesh, has been found to result in significantly improved convergence rates [1 1,13]. The use 
of locally varying smoothing coefficients has the effect of making the scheme more implicit in 
the direction normal to the boundary layer, or normal to the mesh stretching direction, and less 
implicit in the tangential direction. The implementation of implicit residual smoothing with 
locally varying coefficients on unstructured meshes is accomplished by rewritting equation (39) 
as 


Ri = Ri + (40) 

where 1; and q now represent the directions tangential and normal to the local mesh stretching 
vector, as described previously. By analogy with the structured mesh case [11], and making 
use of equation (25), the smoothing coefficients are taken as 


= max 
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CFL* s+1 

where CFL and CFL* are the Courant numbers of the smoothed and unsmoothed schemes 
respectively, s denotes the magnitude of the mesh stretching vector, and <|> is given by equation 
(20). Since equation (40) now contains a directionally scaled Laplacian, it can be discretized on 



-9- 


an unstructured mesh in a manner analogous to that employed for the directionally scaled dissi- 
pation operator, as given in equation (28). For economy, the resulting set of algebraic equa- 
tions are solved only approximately by performing two Jacobi iterations. 

4.3. Multigrid Algorithm 

The idea of a multigrid algorithm is to accelerate the convergence of the fine mesh solu- 
tion by efficiently damping out the low-frequency error components by means of time-stepping 
on coarser meshes. A multigrid method for unstructured meshes has previously been 
developed for inviscid flow calculations [10]. It assumes the various coarse and fine meshes of 
the sequence to be completely independent from one another, and computes the patterns for 
transferring the flow variables, corrections and residuals back and forth between the various 
meshes in a preprocessing operation, where an efficient tree-search algorithm is employed. For 
viscous flow calculations, a full multigrid (FMG) algorithm is employed, where the initial flow 
field on the fine grid is obtained by interpolating a flow solution which has been converged on 
the previous coarser grid with a small number (10 to 20) of multigrid cycles. Better conver- 
gence and additional robustness can also be obtained if the previously employed V-cycle is 
replaced by a W-cycle, which performs one time step on each mesh when proceeding from fine 
to coarse meshes, and no time stepping but merely prolongation of the corrections when 
proceeding from coarse grids to fine grids. It also proves useful to implicitly smooth the 
corrections after the prolongation operation, when proceeding from coarse to fine meshes. The 
constant coefficient implicit smoothing operator of equation (39) is employed for this operation, 
using a value of e=0.2, and the resulting equations are solved approximately using two Jacobi 
iterations. 

5. RESULTS 

The intent of this work is to provide a validation of the basic algorithm described above 
for triangular meshes, and to demonstrate that accurate and efficient solutions can be obtained 
on triangular meshes with highly stretched elements. This is best accomplished by computing 
solutions with the present scheme on triangular meshes which are directly derived from struc- 
tured quadrilateral meshes, and comparing the accuracy and efficiency of these solutions with 
those obtained on equivalent quadrilateral meshes with proven structured-mesh Navier-Stokes 
solvers [11,12,13], 

5.1. Low Reynolds Number Cases 

The first series of test cases involve very low Reynolds number flows over a NACA0012 
airfoil which have been computed by various authors for the GAMM workshop on the Solution 
of Compressible Navier-Stokes Flows [15]. For these cases, the thin-layer assumption does not 
hold, and the flow is dominated by viscous effects, thus providing a means of validating the 
discretization of the full Navier-Stokes viscous terms implemented in this work. The mesh 
employed for these calculations is depicted in Figure 4. It contains 20,800 points and 41,600 
triangles, and is derived from a 320x64 structured quadrilateral C-mesh with 192 points on the 
airfoil, and 64 points in the wake. The far-field boundary is located 15 chords out from the air- 
foil, and the mesh spacing in the normal direction at the wall is 0.002 chords, resulting in rela- 
tively low cell aspect ratios of the order of 10:1 on the airfoil surface, and 100:1 in the wake 
region. For all these cases, a constant temperature wall boundary condition is prescribed along 
the airfoil surface, where the temperature is taken as the adiabatic free-stream temperature. 

In the first test case, the Mach number is 0.8, the incidence is 10°, and the Reynolds 


number is 73. The Mach number contours of the computed solution are depicted in Figure 5, 
where a rapid growth of the boundary layer along the upper airfoil surface is observed, and 
locally supersonic flow is attained only in a small pocket outside the edge of the viscous layer 
on the upper surface. This low Reynolds number flow provides a test of the scheme near the 
Stokes limit of extremely viscous flow. The computed flow field pattern and the lift and drag 
values are within the same range as the results reported in the workshop [15]. A reduction of 
the density residuals of 4 orders of magnitude over 200 multigrid cycles was achieved for this 
case, employing 5 meshes in the multigrid sequence, as shown in Figure 8. 

In the next test case, the Mach number is 0.8, the incidence is 10°, and the Reynolds 
number is increased to 500. The Mach number contours of the computed solution are given in 
Figure 6. A slower boundary layer growth is observed for this higher Reynolds number case, 
accompanied by a stronger leading-edge expansion and an increased region of supersonic flow. 
The recompression from supersonic to subsonic flow appears to occur gradually along the 
upper edge of the viscous layer. Separation occurs on the top surface of the airfoil, and a large 
wake of low-velocity recirculating flow occurs downstream of the airfoil. Nevertheless, rapid 
convergence is achieved with the multigrid algorithm, as shown in Figure 8, where a reduction 
of 10 orders of magnitude of the density residuals is achieved in 200 cycles. The computed 
values of lift and drag for this case compare well with those reported in [11,15]. 

The third case consists of a supersonic low-Reynolds-number flow, where the Mach 
number is 2.0, the incidence is 10°, and the Reynolds number is 106. This represents a standard 
test case which has received wide attention in the literature, and for which experimental data is 
available. The density contours of the computed flow field are depicted in Figure 7, where a 
strong bow shock is observed, which tends to weaken in the far-field due to curvature. These 
computed density contours compare qualitatively with the experimental density contours and 
numerical solutions given in [5,11,15], and the computed lift and drag values are within the 
same range as those reported in these references. For this case, the density residuals were 
reduced by 5 orders of magnitude over 200 multigrid cycles, as shown in Figure 8. 

5.2. Flat Plate Boundary Layer 

An assessment of the accuracy of the scheme may be performed by examining the ability 
of the method to reproduce the well-known compressible boundary-layer solution over a ther- 
mally insulated flat plate. The mesh employed for the boundary-layer calculation is shown in 
Figure 9. It represents a triangulation of a stretched cartesian grid previously employed for 
computing the same problem with a structured mesh solver [16], The mesh contains 24 points 
ahead of the plate, 48 points along the plate in the streamwise direction, and 80 points in the 
normal direction. The upstream boundary is located two plate lengths ahead of the leading 
edge, and the upper far-field boundary is located at a distance of 2.6 plate lengths. The mesh 
points are clustered in the streamwise direction near the leading edge of the plate in order to 
better resolve the stagnation point flow in this region. The mesh point spacing at the wall is 
0.0016 plate lengths, resulting in elements of aspect ratio 50:1 near the trailing edge of the flat 
plate. The computations were performed for a Mach number of 0.8, and a Reynolds number 
based on the plate length of 5000. An exact analytical solution for this flow may be obtained 
by an application of the Howarth-Dorodnitsyn transformation to the incompressible Blasius 
similarity solution [17]. The transformation consists of a rescaling of the coordinate direction 
normal to the plate as a function of the local density variation through the layer: 
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Comparison of computed and exact boundary-layer profiles at the station x=0.6, for a plate of 
length unity, where the Reynolds number based on x is 3000 are shown in Figures 10 through 
12. The normalized streamwise and normal velocities, as well as the shear stress across the 
layer are plotted versus the similarity coordinate r|, which varies from 0 to 1 through the layer, 
and is given by 



where Y is the transformed vertical coordinate given by equation (43). Excellent agreement 
between the computed and exact profiles of streamwise and normal velocity is observed from 
Figures 10 and 11. From Figure 12, good correlation between the computed and exact shear 
stress across the layer is observed. The slight overprediction of the wall shear stress observed 
in this figure, at q - 0, could be systematically reduced by further refinement of the grid. The 
similarity property of the solution was also verified by examining the profiles at various 
different stations along the length of the plate. Good agreement was observed except for sta- 
tions close to the leading edge, where effects of the stagnation point flow are still present, and 
for stations directly adjacent to the outflow boundary. The skin friction along the plate is plot- 
ted in Figure 13, showing good agreement between computed and exact solutions except in the 
aforementioned regions. 

For the present calculations, the k 2 dissipation coefficient was set to zero since the flow is 
subcritical. The value of the K 4 coefficient was taken as 1/256, which resulted in artificial dissi- 
pation terms which were roughly 2 orders of magnitude smaller than the physical viscous terms 
in most regions of the boundary layer. Thus, for these values of the dissipation coefficients, the 
present mesh resolution, which from Figures 10 through 12 can be seen to yield approximately 
20 points in the layer, appears to be sufficient for accurately resolving the boundary layer. A 
reduction of 7 orders of magnitude of the density residuals was achieved over 300 multigrid 
cycles for this case, employing a sequence of 4 meshes, as shown in Figure 14. Similar accu- 
racy could be obtained when the Reynolds number was raised to 50,000, and the mesh spacing 
at the wall was reduced to 0.0005, thus increasing the aspect ratios of the cells, but retaining 
the same number of mesh points in the boundary layer. A somewhat slower convergence rate 
was observed in this case, resulting in a reduction in the residuals of 4 orders of magnitude 
over 300 multigrid cycles, as shown in Figure 14. 

53. Symmetric Laminar Airfoil Case 

The final test case consists of a NACA0012 airfoil at 0 ° incidence, with a freestream 
Mach number of 0.5, and a Reynolds number of 5000. The thermally insulated wall boundary 
condition is applied by prescribing zero heat flux across the airfoil surface. This represents a 
well documented laminar test case which has been computed independently with various struc- 
tured grid codes [11,12,13]. The Reynolds number for this case approaches the upper limit for 
steady laminar flows prior to the onset of turbulence. For this case, separation occurs near the 
trailing edge, and a small symmetric recirculation bubble is formed in the trailing-edge and 
near-wake region. The mesh employed for this case is derived from a 320x64 structured qua- 
drilateral C-mesh with 192 points on die airfoil, and 64 points in the wake, as is depicted in 
Figure 15. It is similar in nature to the mesh of Figure 4, with the exception that increased 
stretching is applied near the airfoil surface and in the wake for better resolution of the thin 
viscous regions. The normal mesh spacing at the wall is 0.0002 chords, resulting in cells with 
aspect ratios of the order of 100:1 along the airfoil. In the wake region, the element aspect 
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ratios were limited to 100:1. Figure 16 depicts the computed Mach number contours in the 
flow field, where the thin boundary-layer and wake regions are visible, and the recirculation 
bubble appears as a region of low Mach number flow. A plot of the velocity vectors near the 
airfoil trailing edge, as shown in Figure 17, clearly shows the recirculatory nature of the flow 
in this region. Plots of surface pressure and skin friction distributions are given in Figures 18 
and 19 respectively. For this subcritical case, the values of the artificial dissipation coefficients 
were taken as k 2 - 0.0, and 1 C 4 - 1/256. This resulted in artificial dissipation terms which were 
roughly 2 orders of magnitude smaller than the physical dissipation terms in the viscous layer 
regions of the flow. The classic trade-off between accuracy and speed of convergence was 
observed for this case by varying the dissipation coefficient k* from 1/256 to 1/64. Table 1 
gives an estimate of the accuracy of the solution as measured by the computed values of pres- 
sure drag, viscous drag, and separation location, versus values produced by various structured- 
quadrilateral-mesh Navier-Stokes solvers. The variation of the multigrid convergence rate with 
the change in the dissipation coefficient is illustrated in Figure 20. The scheme is robust in that 
it converges efficiently over a wide range of K 4 values, and has been found to remain stable for 
values as low as 1/640. While the fastest convergence rate is achieved for a value of K 4 = 1/64, 
the solution accuracy degrades slightly, the separation point having shifted from 81.4% to 
83.4% chord. For all cases, convergence to engineering accuracy could be achieved in less than 
100 multigrid cycles, which requires roughly 7 minutes of CRAY-2 CPU time. Finally, the 
same test case was run on a similar mesh with a normal spacing of 0.00002 chords at the wall, 
with cell aspect ratios of the order of 1000:1 on the airfoil surface and in the wake. For a K 4 
value of 1/64, a reduction of the density residuals of 4 orders of magnitude over 200 multigrid 
cycles was achieved, illustrating the robustness of the code in dealing with the extremely 
stretched elements which are necessary for solving higher Reynolds number turbulent flows. 

6 . CONCLUSIONS 

A new Navier-Stokes solver for use on unstructured triangular meshes has been validated 
by comparing various laminar flow results about simple geometries with well established 
numerical and analytical solutions. The accuracy and convergence efficiency of the present 
scheme were found to be competitive with various well known structured quadrilateral-mesh 
viscous-flow solvers for laminar flow cases. The present code requires approximately 0.19 jc 1(T 3 
sec/node/multigrid cycle of CPU time on a CRAY-2 supercomputer, which represents three to 
four times the computational effort required by equivalent structured-mesh codes. This lower 
computational efficiency is due in large part to the gather-scatter operations required in 
unstructured-mesh algorithms. In future work it will be shown how this solver can be applied 
to arbitrarily complex configurations which are not easily handled by structured-mesh solvers, 
and how the efficiency and accuracy can be improved by the use of adaptive meshing tech- 
niques. A turbulence model for use on unstructured meshes will also be sought for higher Rey- 
nolds number calculations. 
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Figure 1 

Domain of Influence of Node P and Equivalent Control Volume 
for a Finite- Volume Approximation to the Convective Terms 


Figure 2 

Equivalent Control Volume for a Finite-Volume 
Approximation to the Viscous Terms 
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Figure 4 

Stretched Triangular Mesh about a NACA 0012 Airfoil Employed 
for the Low-Reynolds Number Calculations 
Number of Nodes - 20,800, Number of Triangles - 41,600 
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Figure 6 

Mach Number Contours of Computed Solution and Calculated 
Lift and Drag Coefficients due to Pressure Forces for Case 2: 
Mach - 0.8, Re - 500, Incidence - 10® 
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Figure 7 

Density Contours of Computed Solution and Calculated 
Lift and Drag Coefficients due to Pressure Forces for Case 3: 
Mach = 2.0, Re = 106, Incidence - 10" 
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| Figure 8 

Convergence Rate as Measured by the RMS Average of the Density 
Residuals versus the Number of Multigrid Cycles on the Finest Mesh 
for the three Low Reynolds Number Cases: 

Case 1: Mach - 0.8, Re - 73, Incidence * 10° 

| Case 2: Mach * 0.8, Re = 500, Incidence = 10° 

Case 3: Mach - 2.0, Re « 106, Incidence = 10® 
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Figure 9 

Stretched Triangular Mesh Employed for the Flat-Plate Boundary-Layer Calculation 
Number of Nodes - 5913, Number of Triangles - 1 1824 
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Figure 10 

Comparison of Computed and Exact Streamwise Velocity in the Boundary 
Layer at /?£*=3000 in Terms of Similarity Coordinates 
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Figure 13 

Comparison of Computed and Exact Skin Friction Along the Plate 
for a Reynolds Number based on the Plate Length of 5000, and a Mach Number of 0.8 
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Number of Cycles 


Figure 14 

Convergence Rate for Flat Plate Boundary Layer Calculation at Mach - 0.8 
Case 1 : RE L = 5000 Ey walf = 0.0016 

Case 2 : RE L = 50,000 Ay^p 0.0005 
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Figure IS 

Finest Mesh Employed for die NACA0012 Airfoil Calculation at Mach - 03 , Re - 9000 
Number of Nodes - 20,800, Number of Triangles - 41,600 
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Figure 16 

Mach Number Contours of the Computed Solution for Viscous Flow 
Past a NACA 0012 Airfoil, Mach - 0.5, Re - 5000, Incidence - 0 ° 
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Elemt#: 1 Nnode:21440 Ncyc: 200 

cl : 0.0000 cd : 0.0229 cm :-0.0034 


NACA 0012 AIRFOIL (VISCOUS GRID) 

*** PROGRAM NS72 5 STAGE SCHEME *** 

MACH= 0.50 ALPHA= 0.00 RE=0.50E4-04 MGCYC= 2 
VIS0= 1.00 VIS 2 = 0.00 VIS4 = 0.25 


Figure 18 

Computed Surface Pressure Distribution and Calculated Pressure Force 
Coefficients for Flow Past a NACA 0012 Airfoil 
Mach = 0.5, Re = 5000, Incidence * 0° 



Skin Friction Coefficient 
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Elemt#: 1 
cl : 0.0000 


Nnode:21 440 
cd : 0.0332 


Ncyc: 200 
cm : 0.0000 


Figure 19 

Computed Skin Friction Distribution and Calculated Viscous Force 
Coefficients for Flow Past a NACA 0012 Airfoil 
Mach - 0.5, Re - 5000, Incidence - 0* 
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Figure 20 

Convergence Rate as Measured by the RMS Average of the Density Residuals versus the Number 
of Multigrid Cycles on the Fine Mesh for Various Values of the Artificial Dissipation Coefficient K 4 

Mach » 0.5, Re = 5000, Incidence = 0 " 
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NACA 0012 

Mach = 0.5, Re = 5000, 

a = 0 ° 


METHOD 

GRID 

CD p 

CD, 

Separation Point 

Triangle Scheme K* = —77 
256 

320 x64 


0.0332 

81.4% 

Triangle Scheme K 4 = — 7 - 

1 2o 

320 x 64 


0.0336 

82.4% 

Triangle Scheme K 4 = -77 
64 

320 x 64 


0.0344 

83.4% 

Cell-Centered Scheme from [11] 

320 x 64 


0.0337 

81.9% 

Cell- Vertex Scheme from [13] 

256 x 64 


0.0327 

81% 

Cell-Centered Scheme from [13] 

256 x 64 


0.03301 

80.9% 

Cell-Centered Scheme from [13] 

512 x 128 

0.02235 

0.03299 

81.4% 


Table 1 

Comparison of Pressure and Viscous Drag Coefficients and Location of Separation 
Point Computed by the Present Scheme for Various Values of the Dissipation 
Coefficient with Predictions from Various Structured-Quadrilateral-Mesh Codes 
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